Hydrodynamics and rheology of active polar 
filaments 



Tanniemola B. Liverpool 1 and M. Cristina Marchetti 2 

1 Department of Applied Mathematics, University of Leeds, Woodhouse Lane, 
Leeds LS2 9JT, UK t.b.liverpool@leeds.ac.uk 

2 Physics Department, Syracuse University, Syracuse, NY 13244, USA 
mcmOphy . syr . edu 

Summary. The cytoskeleton provides eukaryotic cells with mechanical support and 
helps them perform their biological functions. It is a network of semiflexible polar 
protein filaments and many accessory proteins that bind to these filaments, regulate 
their assembly, link them to organelles and continuously remodel the network. Here 
we review recent theoretical work that aims to describe the cytoskeleton as a polar 
continuum driven out of equilibrium by internal chemical reactions. This work uses 
methods from soft condensed matter physics and has led to the formulation of 
a general framework for the description of the structure and rheology of active 
suspension of polar filaments and molecular motors. 



1 Introduction 

Cells are living soft matter. They are composed of a variety of soft materials, such as 
lipid membranes, polymers and colloidal aggregates, often constrained to a reduced 
spatial dimensionality and geometry. It is then reasonable to expect that the dynam- 
ics and interactions of these constituents that control cell function takes place on the 
same time and energy scales as those of synthetic soft materials. Life adds, however, 
a new feature not found in traditional soft matter: the constant flow of energy and 
information required to keep living organisms alive. This new feature makes cells of 
particular interest to physicists as understanding the behavior of active living matter 
requires the development of new theoretical concepts and experimental techniques. 

The eukaryotic cell cytoskeleton is a perfect example of this novel type of active 
material. The cytoskeleton allows the cell to carry out coordinated and directed 
movements such as cell crawling, muscle contraction, and all the changes in cell 
shape in the developing embryo [1]. The cytoskeleton also supports intra-cellular 
movements such as the transport of organelles in the cytoplasm and segregation of 
chromosomes during cell division [2] . It is highly inhomogeneous, with a large variety 
of different dynamical supramolecular structures. Examples are contractile elements 
like stress fibres, or the contractile ring in mitosis, or astral objects like the mitotic 
spindle which forms during cell division [1, 2]. 
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Self-assembled filamentous protein aggregates play an important role in the me- 
chanics and self-organization of the cytoskeleton. In addition, a number of other pro- 
teins interact with them and modulate their structure and dynamics. Cross-linking 
proteins bind to two or more filaments together to form a dynamical gel. Molecular 
motor proteins bind to filaments and hydrolyze nucleotide Adenosine triphosphate 
(ATP). This process, coupled to a corresponding conformational change of the pro- 
tein, turns stored chemical energy into mechanical work. Capping proteins modulate 
the polymerization and depolymerization of the filaments at their ends. 

A key question is how the elements of the cytoskeleton cooperate to achieve its 
function. To what extent is there a 'cellular' brain and how closely does it control 
cellular mechanisms? How much of a role does spontaneous self-organization driven 
by general physical principles play? 

Much of the recent progress in the understanding of the complex structures 
and processes that control the behavior and function of the cytoskeleton has been 
linked to the development of new biophysical probes allowing an unprecedented view 
of sub-cellular processes at work. Mechanical probes such as optical and magnetic 
tweezers [3], atomic force microscopes [4] and micropipettes probe the response of the 
elements of the cytoskeleton to locally applied forces. Visualization techniques using 
fluorescence microscopy, e.g. fluorescence imaging with one-nanometer accuracy [5] 
or single-molecule high-resolution co-localization [6] based on organic dyes allow 
one to follow the dynamics of single molecules inside living cells, (in-vivo) giving 
insights into the microscopic processes underlying cellular dynamics. Many of these 
experimental developments are reviewed in this book. 

Because of the large number of unknown components, it is also of interest to 
study simplified systems consisting of a smaller number of well characterized ele- 
ments, in-vitro. This has led to a number of experimental biophysical studies of puri- 
fied solutions of cytoskeletal filaments and associated molecular motors that have es- 
tablished that motor-induced activity drives the formation of a variety of spatially in- 
homogeneous patterns, such as bundles, asters and vortices [7, 8, 9, 10, 11, 12, 13, 14]. 
These are reminiscent of some of the supramolecular structures present in the cy- 
toskeleton [15, 16]. The mechanical properties of filament-motor systems has also 
been studied showing qualitative differences from passive filament suspension. Be- 
cause of the controlled nature of their preparation and the detailed knowledge of 
their constituents, in vitro studies are particularly amenable to a quantitative de- 
scription using techniques from theoretical physics. In this review we will be mostly 
concerned with describing the behavior of such simplified systems on large time 
scales (times >microsecond) where the atomistic details are not important and a 
coarse-grained phenomenological description may suffice. 

The reductionist viewpoint typified by this approach also has its drawbacks. A 
simplified system necessarily can provide only a subset of the phenomena observed 
in living cells since only a small fraction of the components are present. A choice 
must also be made of which simplified system to study as different combinations of 
components may give different or similar behavior. This choice must of course be 
heavily influenced by previous experiments [7, 8]. A living cell is a highly optimized 
complex system of interacting agents with the ability to modulate its response to 
complex changes in its environment. This complexity will be missing from simple 
mechanical models described here. There is some hope, however, that this complexity 
can eventually be combined with the physical picture emerging from the approach 
we present here to give a more complete "biophysical" picture of motility in cells in 
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which the laws of physics provide important constraints on the possible "system" 
dynamics. Finally, even within our limited frame of reference, we will also make a 
number of simplifying assumptions in developing the models. Some of the important 
physical phenomena ignored here, such as active polymerization, treadmilling [19, 20] 
and filament flexibility [21, 22, 23, 24, 25, 26, 27], will be incorporated in future work. 

We first review some recent theoretical approaches to describe active filament 
suspensions. We then describe some of our current work and give perspectives for 
the future. 



2 Theoretical modeling of active systems 

There have been a number of recent theoretical studies of the collective dynamics of 
mixtures of rigid filaments and motor clusters. First and most microscopic, numeri- 
cal simulations with detailed modeling of the filament-motor coupling have yielded 
patterns similar to those found in experiments [11, 13], including vortices and asters. 
These simulations modeled the filaments as elastic rods with motor clusters being 
parameterized by three binding parameters, the on and off rates and the off-rate 
at the plus end of the filament. It was found that the rate of motor unbinding at 
the polar end of the filaments plays a crucial role in controlling the vortex to aster 
transitions at high motor densities [12]. 

A second interesting development has been the proposal of 'mesoscopic' mean- 
field kinetic equations first studied in one dimension [28, 29], where the effect of 
motors is incorporated via a motor-induced relative velocity of pairs of filaments, 
with the form of such a velocity inferred from general symmetry considerations. 
Kruse and collaborators [30, 31, 32] proposed a one dimensional model of filament 
dynamics and showed the existence of instabilities from the homogeneous state to 
contractile states [30] and traveling wave solutions [32]. We generalized the kinetic 
model to higher dimensions [33, 34] and used it to classify the nature of the ho- 
mogeneous states and their stability [35, 36]. Related kinetic models have also been 
discussed by other authors [37, 38, 39, 40]. 

Finally, phenomenological continuum theories have been proposed where the 
mixture is described in terms of a few coarse-grained fields whose dynamics is in- 
ferred from symmetry considerations [41, 43, 42, 44, 45, 46, 47, 48, 49, 50, 51]. 

Lee and Kardar [43] proposed a simple hydrodynamic model for the coupled dy- 
namics of a coarse-grained filament orientation and the motor concentration, ignor- 
ing fluctuations in the filament density. These authors argued that filament growth 
by polymerization provides a mechanism for an instability of the system from an 
isotropic to an oriented state [43, 42], with large-scale aster and vortex structures. 
They obtained a phase diagram for the system showing a transition from vortices 
to asters. This model was subsequently generalized by Sankararaman et al. [44] to 
include varying populations of bound and free motors, as well as an additional cou- 
pling of filament orientation to motor gradients. The effects of boundary conditions 
on the steady states of the system was also studied numerically. 

A phenomenological hydrodynamic description for polar gels and suspensions 
including momentum conservation has been discussed by several authors [45, 46, 
47, 48]. These equations generally consider incompressible suspensions and incorpo- 
rate momentum conservation in the Stokes approximation, by assuming the form 
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of constitutive equation for the suspension's stress tensor on the basis of symmetry 
consideration. The coupling of flow and polar order is described via an equation for 
the local polarization of the suspension. This model has been used to identify the 
non-equilibrium defect structures that can occur in the polar state [45, 46] and to 
analyze the behavior of an active polar suspension in specific geometries [52, 47]. 
In particular, it was shown that the interplay between order and activity can yield 
a spontaneous flowing state for a solution near a wall [47]. Closely related hydro- 
dynamic models have been used to describe generically the collective dynamics of 
self-propelled particles in solution, such as swimming bacterial colonies, in both ne- 
matic and polar states [49, 50, 51]. This work builds on earlier work by Toner and 
Tu on hydrodynamic models of flocking, where it was shown that the nonequilib- 
rium nature of internally driven systems allows for novel symmetry breaking phase 
transitions that are forbidden in equilibrium systems with continuous symmetry in 
one and two dimensions [53, 54, 55]. 

The main objective of our work has been to establish the connection between 
microscopic single-polymer dynamics and the phenomenological hydrodynamic mod- 
els by deriving the hydrodynamic equations from a mean-field kinetic equation of 
filament dynamics. In the phenomenological approach the system is described in 
terms of a few coarse-grained fields (conserved densities and broken symmetry vari- 
ables) whose dynamics is inferred from symmetry considerations. The strength of 
this method is its generality. Its drawback is that for systems that are far from 
thermal equilibrium and therefore lack constraints such as those provided by the 
fluctuation-dissipation theorem or the Onsager relations, all the parameters in the 
equations are undetermined. We have bridged the gap between microscopic models 
and continuum theories by deriving the hydrodynamic equations through a sys- 
tematic coarse-graining of the microscopic dynamics. This derivation provides an 
estimate of the various parameters in the equations in terms of experimentally ac- 
cessible quantities. We start with a Smoluchowski equation for filaments in solution, 
where motor proteins are described as active cross-linkers capable of exchanging 
forces and torques between filaments. The active currents arising from such motor- 
mediated exchange of forces and torques are obtained by considering the kinematics 
of two filaments crosslinkcd by a single active protein cluster that can rotate and 
translate at prescribed rates as a rigid object relative to the filaments. The hydrody- 
namic equations are then obtained by suitable coarse-graining of the Smoluchowski 
equation. This method yields a general form of the hydrodynamic equations which 
incorporates all terms allowed by symmetry, yet it provides a connection between 
the coarse-grained and the microscopic dynamics. In a series of earlier publications 
we have described in details the derivation of the hydrodynamic equations for fil- 
aments in a quiescent solvent [33, 34, 56, 35, 36]. Here we generalize this work by 
incorporating the flow of the solvent. This is essential for describing the rheological 
properties of the solution. A brief account of some of the results presented here have 
been given elsewhere [57]. 



3 Hydrodynamics of a solution of polar filaments 

We consider a collection of rigid polar filaments in a viscous solvent. The solution 
forms a quasi-two dimensional film, of thickness much smaller than the length of 
the filaments. Our goal is to study the interplay of order and flow in controlling the 
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phases and the rheology of the system. The filaments diffuse in the solution and can 
be crosslinked by both active and stationary protein clusters. Active crosslinkers 
are small clusters of motor proteins that use chemical fuel as an energy source to 
generate forces and torques on the filaments, sliding and rotating filaments relative 
to each other [2, 12]. In addition, other small proteins, such as a-actinins act as 
stationary crosslinkers and induce filament alignment [1]. 

As in passive solutions of rigid filaments, the large scale dynamics can be de- 
scribed in terms of a set of hydrodynamic equations for continuum fields that relax 
on time scales much longer than microscopic ones. These include the conserved vari- 
ables of the systems, as well as any field associated with broken symmetries. Various 
forms of these equations have been written down phenomenologically by other au- 
thors. What distinguishes our work from these phenomenological approaches is that 
we derive the hydrodynamic equations from a mesoscopic model of coupled motor- 
filament dynamics. This allows us to estimate the various parameters in the hydro- 
dynamic equations (which are undetermined in the phenomenological approach) and 
relate them to quantities that can be controlled in experiments. To make contact 
with the existing literature, we first present the equations and then discuss their 
derivation via coarse-graining of a Smoluchowski equation for rigid rods in a viscous 
solvent. 

The conserved densities in a suspension of interacting filaments (rods) in a 
solvent are the mass densities of filaments (rods) p r (r,t) — mc(r,t) and solvent 
p s (r,t), and the total momentum density g(r, t) = p(r, t)v(r, t) of the solution 
(rods+solvent), with v(r,t) the flow velocity and p(r,t) = p s + p r the total density. 
Here c(r, t) is the number density of rods and m the mass of a rod. The conserved 
densities satisfy conservation laws, given by 

dtp = -V • g , (1) 
d t c = — V • J , (2) 

d t gi + d j (g l g j /p) = djOi, + pF, ext , (3) 

where J(r, t) is the current density of rods and F ext the external force on the sus- 
pension. The stress tensor Oij is the i-th component of the force exerted by the 
surrounding fluid on a unit area perpendicular to the j-th direction of a volume 
element of solution. It includes all forces on a volume of suspension exerted by the 
surrounding fluid. It can be written as the sum of solvent and filament contributions 
as 

°ij = a lj + °lj ■ (4) 
The solvent contribution has the usual form appropriate for a viscous fluid, 

alj = 2?7oWij + (Vb - rio)5xjU k k - SijII s (p) , (5) 

where r/o and rjb are the shear and bulk viscosity of the solvent, 17 a (p) is the pressure 
of the solvent, and Uij is the symmetrized rate of strain tensor, 

U H = \ (dm + d 3 v i) ■ ( 6 ) 

In the low Reynolds number limit we can ignore the inertial terms on the left hand 
side of the solvent momentum equation, Eq. (3). 
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In a solution of long filaments states with liquid crystalline order are possible. 
Polar rods can form polarized and nematic states, both characterized by orienta- 
tional order, but with different symmetries for the order parameters. Polar order 
in a fluid of iV rods is characterized by a vector order parameter or polarization, 
P(r, t), defined by 

JV 

P(r,t) = ^y(2fla(t)*(r-ra(t)), (7) 

a — l 

where v a is the position of the center of mass of the a-th rod and u a is a unit vector 
directed along the polar direction. The angular bracket denote an ensemble average. 
The polarization vector P can be written as 

P(r,i) = P(r,t)p(r,i) , (8) 

where the magnitude of the polarization P(r, t) is the scalar order parameter and 
the unit vector p(r, t) identifies the direction of broken symmetry in the ordered 
state. 

Nematic order is described by the conventional nematic order parameter or align- 
ment tensor, defined as 

JV 

Q«(r,t) = ( Yl (««•«(*)«««(*) - - •■«(«)) • ( 9 ) 

a — l 

The subtracted part ensures that the order parameter vanishes in the isotropic state 
in d dimensions. The alignment tensor Qij is thus a traceless and symmetric second- 
order tensor field, with two independent degrees of freedom in d = 2. For uniaxial 
nematics the alignment tensor takes the form 



Qij = S(r,t) 



n t {r,t) nj (r,t) - ^Sij , (.1.0) 



where S(r, t) is the scalar nematic amplitude and n(r, t) is the familiar nematic 
director. The nematic state has orientational order (S ^ 0) and it is invariant under 
inversion of the director, i.e., for n — * — n. The polarized state, in contrast, is not 
invariant for p — > —p. In a polarized state the alignment tensor Qij is slaved to the 
polarization and acquires a nonzero value, with n = p. 

The dynamical equations for polarization and alignment tensor have the form 
(for simplicity we give the form for d = 2 only) 

A(cP) =cPi(#e,P) -djJij -P< , (11) 
D t (cQij) = cFij (k, Q) -dkJijk - Rij , (12) 

where D t — d t + v • V, k is the rate of strain tensor, Kij = djVi, and 

5 

Fi(n,P) = -LOijPj + XpUijPj - -UkkPi , or 

F(k, P) = i(V x v) x P + A P [Vv + (Vv) T ] • P - j(V • v)P , (13) 
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Fij(K, Q) — — (ulikQkj + ^>jkQki) + -^yUikQkj +UjkQki ~ SijUklQklJ 

--UkkQij + ^(uij - ^SijUkk^ • (14) 

Here Ap and A are the flow alignment parameters in the polarized and nematic 
states, respectively, and 

Uij = ^(diVj - djVi) . (15) 

The low density derivation based on the Smoluchowski equation described below 
gives Ap = 1/2 and A = l/(2£ ), with So the magnitude of the nematic order 
parameter. Since typically in the nematic state even quite far from the I-N tran- 
sition, 5*0 <C 1, we expect A > 1, as required for flow alignment. In addition it is 
well known that deep in the nematic phase, higher order correlations can further 
increase the value of A. In the following we will treat both Ap and A as unknown 
parameters. The first terms on the right hand side of Eqs. (11) and (12) general- 
ize the convective derivative on the left hand side of the equation to the case of 
long, thin rods. These are standard terms in nematohydrodynamics that have been 
derived from a microscopic model before [58]. Different values for some of the nu- 
merical coefficients are reported in the literature, depending on the closure scheme 
used in evaluating various angular averages. The second and third terms on the 
right hand side of Eqs. (11) and (12) represent translational and rotational currents, 
including contributions from diffusion, excluded volume, and both stationary and 
active cross- linkers. The relaxation of the order parameters is controlled by the rota- 
tional currents and is non- hydro dynamic. In contrast, the relaxation of the broken 
symmetry variables p and n is controlled by hydrodynamic Goldstone modes, as 
appropriate in ordered states. 

The long wavelength description of the solution is then given by the five equa- 
tions (1-3) and (11-12). To close the hydrodynamic equations we must derive the 
constitutive equations for the fluxes (J, JV,, Jijk), the rotational currents (Ri and 
Rij), and the filament contribution to the stress tensor, <j\j, as functions of the sys- 
tem's properties (density, filament concentration and order parameters) and of the 
driving forces (applied mechanical stresses and activity, as measured by the ATP 
consumption rate). This derivation is carried out below by adapting methods from 
polymer physics appropriate for a dilute solution of rigid rods. Although the specific 
expressions obtained by this method for the parameters in the hydrodynamic equa- 
tions only apply at low concentration of filaments, the structure of the equations is 
general and remains the same at high density. 



4 Derivation of hydrodynamic constitutive equations 

Our goal is to derive the constitutive equations for the various hydrodynamic cur- 
rents and stresses starting from a semi-microscopic model of the dynamics of single 
filaments coupled pairwise by active and stationary crosslinkers. The filaments are 
modeled as rigid rods of fixed length I and diameter b <C I immersed in a viscous sol- 
vent. They diffuse independently in the solvent and interact via excluded volume. In 
addition, filaments can be coupled pairwise by both stationary and active crosslink- 
ers that generate additional active currents. Active crosslinkers are described as rigid 
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links that can walk along the filaments towards the polar end at a prescribed rate 
u(s) proportional to the rate of ATP consumption. Generally u(s) varies with the 
point s of attachment along the filament (0 < s < I). Both active and stationary 
cross-links also mediate the exchange of torques between the filaments by acting as 
torsional springs of prescribed stiffness, k. Our goal is to obtain a coarse-grained de- 
scription of the system where all the parameters in the hydrodynamic equations are 
characterized in terms of u(s), n, and the density of crosslinkers. Collective effects 
arising from multiple crosslinkers are neglected and the density of crosslinkers is 
assumed constant for simplicity. We also neglect the dynamics of crosslinkers bind- 
ing and unbinding which occurs on faster time scales than those of interest here, so 
that we can treat a constant fraction of them as bound. The dynamics of crosslinkers 
binding and unbinding was considered for instance in Ref. [44] and it was found that 
varying the rates of motor binding and unbinding did not affect the nonequilibrium 
steady states of the active solution. The derivation of the active contributions to the 
various fluxes has been presented elsewhere [36] and will be summarized here for 
completeness. We also present novel results on the evaluation of the filament contri- 
bution to the stress tensor up to terms of first order in gradients of the hydrodynamic 
fields. 

To proceed, we also make a series of simplifying assumptions on the dynamics 
of the solution. First, we assume that the friction between filaments and solvent is 
large and the filaments move at the flow velocity v = g/p of the solution. In many 
fluid mixtures internal friction mechanisms are so strong that the flow velocities of 
the two components relax on microscopic time scales to the common value v. There 
are situations, however, where the relaxation time of the relative momenta of the 
two species is slow enough to have a significant influence even on hydrodynamic time 
scales. In this case a two-fluid description is appropriate and useful. Such a "two-fluid 
model" of the system (rods and fluid background) will be described elsewhere, where 
we will show under which conditions one approaches the one-fluid model (which is 
always the true hydrodynamic limit). 

We also limit ourselves to the case of incompressible solutions, with p — p s + p r = 
constant, which requires 

V-v = 0. (16) 

Finally, we neglect fluid inertial effect compared to the frictional forces between the 
colloidal rods and the solvent. In this limit the momentum equation (3) reduces to 
the Stokes equation 

= - P Fr , (17) 

or, in the absence of external forces, 

rjoV 2 Vi - dill 3 = —djCJij . (18) 

Equation (18) shows that the flow velocity of the solution is determined by the 
stress introduced by the filaments. In turn, the forces that the filaments exert on 
each other and on the solvent depend on the flow of the suspension in which they 
are immersed and the problem must be solved self-consistently. 

The dynamics of a dilute suspension of rods in the presence of a macroscopic 
flow field v(r) can de described by the Smoluchowski equation for the probability 
distribution c(r, u, t) of rods with center of mass at r and orientation u at time t . The 
Smoluchowski equation describes the mean-field Brownian dynamics of extended 
colloidal particles at low Reynolds number, under the assumption that the particles 
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velocities have equilibrated on microscopic time scales to a local Maxwell-Boltzmann 
distribution at a temperature T a [59, 60]. The effective temperature T a incorporates 
nonthermal noise sources as may arise from fluctuations in motor concentration and 
ATP consumption rate. The Smoluchowski equation is given by 

dtc + V -Jc + K- J c = , (19) 

where R = u x 9 a is the rotation operator. The translational probability current, 
J c , and the rotational probability current, J c , are given by 

Jd = &v t ~ D lJ V J c - T^rc VjUe* + j£ , (20) 

J ci = CU)i - DrTZiC - Dr clliUcx + Jci , (21) 

where w» = eijkUjUidiVk- Also Dij = D^UiUj + D±(Sij — UiUj^j is the translational 
diffusion tensor and D r is the rotational diffusion rate. For a low-density solution 
of long, thin rods D± = Dy/2 = D/2, where D = k B T a ln(l/b)/(2ivri l), and D r = 
QD/l 2 . The potential U cx incorporates excluded volume effects which give rise to 
the nematic transition in a solution of hard rods, ft can be written by generalizing 
the Onsager interaction to inhomogeneous systems as fcsT a times the probability of 
finding another rod within the interaction area of a given rod. In two dimensions 
this gives 

E/ex(ri,ui) =k B T a I du 2 I |ui x u 2 | c(n +£,u 2 ,t) , (22) 



where Sj, with — 1/2 < s, < 1/2, parameterizes the position along the length of 
the i-th filament, for i = 1,2, and ... = J*^J 2 dsi... = (..) Si - The filaments are 
constrained to be within each other's interaction volume, i.e., in the thin rod limit 
b <ti I considered here, have a point of contact. The factor |ui x U2I represents 
the excluded area of two thin filaments of orientation ui and u.2 touching at one 
point [61]. Finally, £ = r-2 — ri ~ uisi — U2S2, is the separation of the centers of 
mass of the two rods. The translational and rotational active current of filaments 
with center of mass at ri and orientation along ui are written as 



J*(ri,ui) = c(n, ui, t)b 2 m / / v a (l; 2)c(ri + £, u 2 , 

J U2 J s l s 2 



t) , (23) 



Jc{n,ui) =c(n,ui,t)6 2 m / / Wo (l; 2)c(n + ^, u 2 , t) , , (24) 

where m is the density of bound crosslinkers and (1; 2) = (si, ui; s 2 , u 2 )- Finally, 
v a (l;2) and oj a (l;2) are the translational and rotational velocities, respectively, 
that filament 1 acquires due to the crosslinker-mediated interaction with filament 2, 
when the centers of mass of the two filaments are separated by £ (see Fig. 1). 

The derivation of the form of the active velocities in terms of motor parameters 
(the stepping rate u(s) and the torsional stiffness ft) has been discussed in detail 
elsewhere [56, 36]. The angular velocity is 

u« = 2 [7P + 7atp (ui ■ U2)] (ui x u 2 ) , (25) 
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Fig. 1. The geometry of overlap between two interacting filaments of length I cross- 
linked by an active cluster. The cross-link is a distance si, (S2) from the center of 
mass of filament 1(2). The distance between centers is £ = V2 — ri = sini — S2n2. 

with 7p and 7jvp motor-induced rotation rates due to polar and nonpolar crosslink- 
ers, respectively (see Fig. 2) . The motor-induced translational velocity has the form 

Va(l; 2) = |v r + Vm, with [36] 




(a) (b) 

Fig. 2. Polar and nonpolar clusters interacting with polar filaments. Assuming that 
clusters always bind to the smallest angle, polar clusters bind only to filaments in 
configuration (a) while non-polar clusters bind to both configurations equally. 



Vr = -(U 2 

V m = A (u 2 + Ul) + B (u 2 - Ul) , 

where a, = a(l + ui ■ Us) and j3 = (3(1 + ui ■ u 2 ). The expressions for A and B 
have been obtained in [36] using momentum conservation. For long thin rods with 
C_l = 2£|| = 2^, to leading order in ui • U2, we find A = — [/3 — a(si + S2)/2]/12 and 
B = a(si - s 2 )/24. 
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The rotational rates, 7p and 7atp, and the active velocities a and (3 can be 
related to the torsional stiffness n of the crosslinkers and to the rate u(s) at which 
a motor cluster attached at position s steps along a filament towards the polar end. 
This rate will in general depend on the point of attachment s, due for instance to 
crowding or stalling of motors near the polar end. The mean (averaged along the 
filament) stepping rate uo = { u ( s ))/l is simply proportional to the mean rate Ratp 
of ATP consumption via a the characteristic step length, which we take of order 
of the thickness b of the filaments, uo ~ W?atp- We emphasize, however, that in 
general the stepping rate u(s) (and other active parameters) may be a non-linear 
and possibly even non-monotonic function of the rate of ATP consumption, Ratp- 

In our model there are three coupled mechanisms for crosslinker-induced filament 
dynamics, described by the parameters a, f3 and the rotational rates, 7p and 7jvp- 
The first is the bundling of filament of similar polarity at a rate a given by [56] 

a = / ' ^ 7 < s ) « Mb/l) , (26) 

J -1/2 1 1 

where the last approximate equality applies in situations where u(s) exhibits strong 
spatial variations on length scales of order b, as may arise for instance from motor 
stalling at the polar end [56]. It is apparent from Eq. (26) that a = if u(s) is con- 
stant. Bundling is driven by the contractile nature of motor clusters and in our mean 
field model requires spatial inhomogeneities in the rate at which motors step along 
the filaments. As we will see below, it tends to build up density inhomogeneities and 
is the main pattern-forming mechanism. The second mechanism of motor-induced 
dynamics will be referred to as "polarization sorting", although in general it involves 
coupled filament rotation and spatial separation of filaments of different polarity. It 
occurs at the rate 

f l/2 ds 

H = / -u(s) = uo , (27) 

J-l/2 1 

and vanishes for aligned filaments. This mechanisms is especially important in the 
polar state where it allows for terms in the hydrodynamic equations corresponding 
to convection of filaments along the direction of mean polarization and it provides 
the mechanism for the transition to a state with spontaneous flow [47]. Finally, 
motor- induced filament rotations occur at rates 7p and 7jvp for crosslinkers that 
preferentially bind to filaments of the same orientation (jp) or regardless of their 
orientation (7jvp). As discussed in Ref. [36], we estimate 

7P ~ 7JVP ~ , (28) 

with £ r = kp,T a ID r a rotational friction. In general both active and stationary 
crosslinkers may induce rotation and be either polar or apolar in nature. In the 
following we will restrict ourselves for simplicity to the case where all polar cross- 
linkers are active motor clusters (of density m a ), while all apolar crosslinkers are 
stationary (of density m s ). In practice we do expect this to be often the case. The 
rotational rate 7p will then depend on ATP consumption, with -yp ~ Ratp, while 
we expect 7jvp to be essentially independent of it. As mentioned above, the various 
active parameters may be non-linear and even non-monotonic functions of Ratp- 
However, these effects will not be considered here. 
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From the Smoluchowski equation, Eq. (19), we obtain the hydrodynamic equa- 
tions for filament concentration, polarization and alignment tensor by truncating 
the exact moment expansion of c(r, u, t) as 

c(r, u, t) = { 1 + 2P(r, t) ■ u + 4Q lJ (r, t)Q t] (u) + . . . } , (29) 

with Qij(u) — UiUj — ^Sij and keeping only the first three moments, 

J du c(r, u, t) — c(r,t) (density), 

J du u c(r, u, t) = c(r, i)P(r, t) (polarization), (30) 
J du Qij (u) c(r, u, t) = c(r, t)Qij (r, t) (nematic order) . 

The details of the calculation, which involves using a small gradient expansion for 
the filament probability distribution and evaluating angular averages, are given in 
[36], where the full expression for the various fluxes and rotational currents are also 
displayed. The resulting hydrodynamic equations in the isotropic and ordered phases 
will be given below. 



5 Stress tensor of an active solution 

In this section we derive the constitutive equation for the filament contribution to 
the stress tensor of an active suspension of polar filaments. An important difference 
as compared to passive solutions is that in active systems stresses can be induced not 
just by externally applied mechanical deformations (yielding k„ / 0), but also by 
motor activity which maintains the system out of equilibrium by supplying energy 
at a rate Ratp- 

In the limit where inertial effects may be ignored (low Reynolds number) and in 
the absence of external forces, momentum conservation is described by Eq. (18), with 
V ■ v = 0. Using standard methods from polymer physics, the filament contribution 
to the divergence of the stress tensor of a dilute suspension of hard, thin rods can 
be written as 

V-«7 r = - JJ^{r-t,u,t)(5((--su)F h {sj) s , (31) 

where T h (s) is the hydrodynamic force per unit length exerted by the suspension on 
a rod at position s along the rod. It arises from interactions with other filaments and 
proteins and with solvent molecules. It depends implicitly upon direct interactions 
between the rods, as well as on hydrodynamic interactions mediated by the solvent. 

The hydrodynamic force density on a rigid rod suspended in a viscous solvent 
can be expressed in terms of the force and torque at its center of mass. A sketch of 
the derivation is given in Appendix A. Further details of similar calculations can be 
found in [61, 62]. We find that the stress due to the filaments can be written in the 
form (to C(V 2 )) 
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V-cr r (r,t) = / c(r,u,t)F h (r,u,t) 

J U 

J u 

In the absence of inert ial effects, the total hydrodynamic force, F h (r,u,t), exerted 
by the suspension on the center of mass of a rod can be found from the condition 
that all forces acting on the rod must balance. The solvent flow field on a given 
segment of a rod is calculated using a decoupling approximation where the hydro- 
dynamic coupling to other segments of the same rod are treated explicitly within 
the Oseen approximation, while the hydrodynamic effects of other rods enter in the 
determination of a self-consistent value for the flow velocity of the solvent, yielding, 

F h (r, u, t) = k B T a V In c + Vt/ ex - F a , (33) 

where — fcsT a Vlnc is the Brownian force, — V[/ C x is the force due to the direct 
interaction of the rod with other rods (in this case, via excluded volume) and F a is 
the active force that can be written as 

F ai = fc(u)J^/c . (34) 

The rod friction tensor dj(u) is proportional to the inverse of the rod diffusion 
tensor Dij(u), with 

&(u) =fe B T a [D- 1 (u)] ij 

= C\\UiUj + C±(6ij - Uiiij) , (35) 

with £|| = 2-Kr/ol/ \n(l/b) and (± = 2£y. Similarly the total hydrodynamic torque is 
given by 

T h (r,u,t) = [k B T a Ti\nc + nU x - t ] x u 

-|uu(u-V).v(r), (36) 

with T a = (r j7c A /c the active torque. The last term on the right hand side of Eq. (36) 
is a viscous contribution to the stress proportional to the velocity gradient. 

The rod contribution to the stress tensor can now be evaluated explicitly using 
the truncated moment expansion for c(r, u, i) given in Eq. (29). When evaluating 
the active contributions to the stress tensor, only terms up to first order in ui ■ U2 
are retained in the active force C(ui) • v a (l;2) exerted by a motor cluster on the 
filament. This approximation only affects the numerical values of the coefficients in 
the stress tensor, not its general form. 

For simplicity, we consider solutions in the presence of a constant velocity gra- 
dient, Kij, and with a uniform mean rate of ATP consumption. We allow for spatial 
inhomogeneities in the filament concentration and orientational order parameters 
and evaluate the stress tensor up to first order in gradients of these hydrodynamic 
fields. The deviatoric part dij = ov, — (l/2)5ijakk of the stress tensor of the filaments 
is 

flr&(r,t) = ff3(r,t) + *y(r,t) , (37) 

with 
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atj = 2k B T a c(l - —)Qij - kBTa — ip.Pj - \ 

\ ClN ! Cip \ I 

, ,2 k B TJ A 2 /4 1 2 \ 

a a 72D V3 j 3 ~ 2 3 J 



+m a b (3 



2 a kBT a l 4 



216D 



^Pi-^V-P - - (i),J>, - <),!■', 



(38) 



where cip = D r /(m a b 2 -ypl 2 ), and cin = cjv/[1 + C]v/ 2 m s 6 2 7jvp/ (4D r )] are the den- 
sities for the isotropic-polarized (IP) and isotropic-nematic (IN) transition, respec- 
tively, at finite density of active polar motor clusters (m a ) and stationary nonpolar 
crosslinkers (m s )[36]. Finally, cjv = 37r/(2/ 2 ) is the density of the IN transition in 
passive systems. There are three types of contributions to the active part of the stress 
tensor. The first consists of the first two terms on the right hand side of Eq. (38). 
These are equilibrium-like terms, in the sense that they have the same structure one 
would obtain in a nematic and polar passive fluid, respectively, with the transition 
densities replaced by their active values. In particular, the first term on the right 
hand side of Eq. (38) should be compared to the corresponding contribution for 

isotropic (c < cn) passive solutions, afj = 2fcsTc^l — -^jQij. The third term is 

a homogeneous nonequilibrium contribution that remains nonzero even for mj = 0. 
This "spontaneous stress" arises from activity and is proportional to the ATP con- 
sumption rate that acts as an additional driving force and can build up stresses even 
in the absence of mechanical deformations. This term is generated by motor-induced 
filament bundling and it is proportional to the bundling rate, a. It would therefore 
vanish in the absence of spatial inhomogeneities in the motor stepping rate. Finally, 
the fourth term contains active contributions proportional to gradients of the po- 
larization (we have omitted here terms of linear order in the gradients containing 
both polarization and alignment tensor. The full expression for the stress tensor can 
be found in Appendix A). These stresses are generated by motor- induced filament 
sorting and are proportional to /3. They are important only in the polarized phase, 
where we expect they will play an important role in enhancing the relaxation of 
longitudinal fluctuations of the filaments and the corresponding relaxation of shear 
via reptation. 

Finally, the viscous contribution to the stress is 



Hi 

48 L2^' J 2 J > ' :i 
+ - {uikQkj + UjkQkij 



(39) 



6 Homogeneous states of a quiescent solution 

We first examine the case of a quiescent suspension, with v = 0. We consider a 
system with a concentration m a of active, polar motor clusters and a concentration 
m s of stationary nonpolar crosslinkers. For convenience we define a dimensionless 
parameter fi a measuring activity as 

fi a = m a b 2 jL ~ m a R A TP , (40) 
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where D r is the rods' rotational diffusion constant and we have assumed that 7p ~ 
Ratp- We also introduce a dimensionless parameter fj, s measuring the effect of 
stationary crosslinkers as 

[is = msb , (41) 

and assume that fi s is essentially independent of the ATP consumption rate. The 
bulk states of the system are determined by the solution of the homogeneous hy- 
drodynamic equations containing only those terms that are of zeroth order in the 
gradients. This is the most coarse-grained description of the system. More detailed 
descriptions that incorporate slowly varying spatial variations can then be devel- 
oped by including gradient terms in the hydrodynamics. The possible homogeneous 
states of the system are obtained as the stationary solution of the homogeneous hy- 
drodynamic equations for filament concentration, polarization and alignment tensor, 
setting all gradient terms equal to zero. In this case the filament concentration is 
constant, c = Co, and only contributions from rotational currents survive in equation 
for the polarization and the alignment tensor, which are given by 



dtPr = -D r [l - flaCo] Pi + D r 



4c /cjv + (p s - 2fj, a )c 



QijPj , (42) 



dtQij = -D f 



4(1 - cq/cn) - Ms c 



(lj +2D r ,i ll cu(P i P j -^SijP 2 ^ , (43) 



where all filament densities are measured in units of I 2 , and cjv = 3tt/2. 

There are three possible homogeneous stationary states for the system, obtained 
by solving Eqs. (42) and (43) with d t Pi = and dtQij = 0. These are: 

isotropic state (I) : Pi = Qij = , 
nematic state (N) : Pi = Qij ^ , 
polarized state (P) : Pi ^ Qij / . 

At low density the only solution is Pi = and Qij — and the system is isotropic 
(I). The homogeneous isotropic state can become unstable at high filament and/or 
motor density, as described below. 

To discuss the instabilities it is convenient to measure time in units of D^ 1 and 
rewrite Eqs. (42) and (43) in a more compact form as 

dtP* = -aiPi + biCoQijPj , (44) 

dtQij = -a 2 Q 1] + 6 2 co (PiPj - ^P 2 ^ . (45) 
The coefficients oi, fei, a 2 , and b 2 are given by 

ai = 1 - fi a co , (46) 

o 2 = 4[l - c /c N - n s co/4] , (47) 

bi = 4/cjv + Us - 2fi a , (48) 

b 2 = 2fi a . (49) 



In the absence of crosslinkers (/x s = /j, a = 0), no homogeneous polarized state 
with a nonzero mean value of P is obtained. There is, however, a transition at the 
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density cjv = 37r/2 from an isotropic state with Qij = for c < cat to a nematic 
state with Qij = So(nirij — \8ij), with n a unit vector along the direction of broken 
symmetry, for Co > cat. The transition is identified with the change in sign of the 
coefficient a 2 of Qij on the right hand side of Eq. (45) . A negative value of a 2 that 
controls the decay rate of Qij signals an instability of the isotropic homogeneous 
state. A mean-field description of such a transition, which is continuous in 2d (but 
first order in 3d), requires that one incorporates cubic terms in Qij in the equation 
for the alignment tensor. Adding a term —aACoQklQklQij to Eq. (45) we obtain 
So = ^ ^-2a 2 /a 4 = ± \/-8(l - c /c N )/a 4 . 

If fi a = 0, but fi s 7^ 0, there is again no stable polarized state. The presence of 
a concentration of nonpolar crosslinkers does, however, renormalize the isotropic- 
nematic (IN) transition, which occurs at a lower filament density given by 

cin = - ■ CN . (50) 
I+HsCn/4 

The presence on nonpolar crosslinks favors filament alignment and shifts cim down- 
ward. It should be noted that this occurs even with a higher effective temperature 
T a . A qualitatively similar result has been obtained in numerical simulation of a two- 
dimensional system of rigid filaments interacting with motor proteins grafted to a 
substrate [64]. The amount of nematic order So is also enhanced by the crosslinkers, 

with So = y / -2a 2 /a 4 cl = ^J-^t( c o/cin - 1). 

If fi a is finite, the system can order in both polarized and nematic homogeneous 
states. The homogeneous isotropic state can become unstable in two ways. As in the 
case ^ = 0, a change in sign of the coefficient a 2 signals the transition to a nematic 
(N) state at the density c/jv given in Eq. (50). In addition, the isotropic state can 
become linearly unstable via the growth of polarization fluctuations in any arbitrary 
direction. This occurs above a second critical filament density, 

cip = — , (51) 

Ha 

defined by the change in sign of the coefficient 01 controlling the decay of polarization 
fluctuations in Eq. (44). For Co > cip the homogeneous state is polarized (P), 
with P / 0. The alignment tensor also has a nonzero mean value in the polarized 
state as it is slaved to the polarization. The location of the boundaries between the 
various homogeneous states is controlled by the relative strength and concentration 
of active polar motor clusters to stationary nonpolar crosslinkers. In order to simplify 
the discussion we fix the value of fis that determined the density of the nematic- 
isotropic transition to /x s = 0, so that the isotropic-nematic transition takes place 
at the density cjv of a suspension of rods with no crosslinkers. One can identify two 
regimes. 

I) If cip > cat, which corresponds to Ha < 1/cjv, a region of nematic phase exists 
between the isotropic and the polar state. At sufficiently high filament and motor 
densities, the nematic state becomes unstable. To see this, we linearize Eqs. (44) and 
(45) by letting Qij — Q°j +SQtj and 8 Pi — Pi. Fluctuations in the alignment tensor 
are uniformly stable for a 2 < 0, but polarization fluctuations along the direction of 
broken symmetry become unstable for 01 < cofciSo/2, i.e., above a critical density 



1 

CArp = — 

Ha 



(52) 
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where R — cisr/cip. The polarized state at Co > cnp has 

P° = Pop. , (53) 
Qij = S P (piPj - Sij/2) , (54) 

with p a unit vector in the direction of broken symmetry and 

C50162 L \oic 5o/ J 
V 2aia 2 

II) For fi a > 1/cjv, c/p > cjv and the polarity of motor clusters renders the 
nematic state unstable at all densities and the system goes directly from the I to the 
P state at c/p, without an intervening N state. The phase diagram has the topology 
shown in Fig. 3. At the onset of the polarized state the alignment tensor is again 
slaved to the polarization field, Qij = (P%Pj ~ \&ijP 2 ) 1 an d P = -Pop. The 

value of Pq is determined by cubic terms in Eq. (44) not included here. 



a 1 



(56) 




For a fixed, but nonzero value of /x s , the phase diagram has the same topology as 
shown in Fig. 3, but with cjv replaced by c/jv given in Eq. (50). The value of \x a where 
the three phases coexist is shifted to a larger value, given by /i a = (1 + /i s cjv/4)/cjv. 

Estimates of the various parameters can be obtained using a microscopic model 
of the motor-filament interaction of the type described in Appendix A. Using pa- 
rameter values appropriate for kinesin (k ~ 10 -22 nm/rad [65]) we estimate jp ~ 
jnp ~ n/Cr = kD t I '(fcsT,,) ~ 10 _1 sec~ 1 , where we used the value D r ~ 10 _2 sec _1 
appropriate for long thin rods in an aqueous solution [66] and T a ~ 300 K. Using 
m a = m a = m, 7p = jnp, I ~ 10/im, b ~ 10 nm, the value of m above which no 
nematic state exist is found to correspond to a three-dimensional crosslinker density 
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of about 0.5 — lf/M and a sample thickness of order lfim. This value is of order of 
the motor densities used in experiments on purified microtubule-kinesin mixtures 
such as those of Ref. [17] , suggesting that the filament solution in this experiments is 
always in the region where the present mean field model predicts a uniform polarized 
state. 

On the other hand, in vitro experiments generally fail to observe states with uni- 
form polarization and report the formation of complex spatial structures. This can 
be understood in the context of the hydrodynamic theory described here by exam- 
ining the dynamics of spatially varying fluctuations of the hydrodynamic fields from 
their uniform value in each state. It has been shown elsewhere that such fluctuations 
become unstable in a wide range of parameters. In both isotropic and ordered states 
the instability arises from filament bundling (controlled by the rate a) that tends 
to build up density inhomogeneities, eventually overtaking diffusion and driving the 
formation of spatially inhomogeneous structures. This instability is described in the 
next section for the isotropic state. The instability of the nematic and polarized 
state is driven by the same physical mechanisms, although the details are more sub- 
tle as in this case one must consider the coupled dynamics of fluctuations in the 
concentration and in the orientational order parameter. A complete description can 
be found in Ref. [36]. 



7 Hydrodynamics of flowing active suspensions 

In this section we display the explicit form of the hydrodynamic equations for a 
suspension of active rods obtained by coarse- graining the Smoluchowski equation, 
as outlined in Sections 4 and 5. Phenomenological forms of these equations have 
already been used by other authors to study the interplay of order and flow in 
active systems in specific geometries [47, 52, 51]. Our work provides a derivation 
of the continuum theory starting from the dynamics of single filaments coupled 
by active crosslinkers and an estimate of the various parameters in the equations 
in terms of experimentally accessible quantities. As discussed in Section 4, we limit 
ourselves to the case of an incompressible suspension and neglect inertial fluid effects. 
In this case the flow velocity v of the suspension is determined from the solution of 
Stokes' equation, 

7?oV 2 v- V77(c,P,Q;k, M ) = -V • * r (c, P, Q; «, ^); , (57) 
with the incompressibility condition 

V • v = . (58) 

The pressure 77 is the sum of solvent and filament contributions, 77 = n s (p) + 77 r , 
and we have introduced the deviatoric stress tensor defined by subtracting out the 
hydrostatic pressure, 77 r = (l/d)o r kk , as 

dij = o~lj - <5 i3 77 r . (59) 

Both isotropic and ordered (polarized and nematic) suspensions will be considered. 
The orientational order of the suspension affects the flow through the dependence of 
the pressure 77 and the rods' contributions to the stress tensor <r r on polarization and 
alignment tensor. The derivation of the constitutive equations for these quantities 
was described in Section 5. 
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7.1 Isotropic state 

In an isotropic suspension the only hydrodynamic variable describing the filaments 
is the concentration, c. Its dynamics is governed by a nonlinear convection-diffusion 
equation 

<9tc + V ■ (vc) = V ■ P(c)Vc , (60) 

where X>(c) is an effective (concentration-dependent) diffusion coefficient, softened 
by active processes. It is given by 

3D 

T>(c) = — (1 +voc) — arh a c , (61) 

with fh a = m a b 2 . The first term on the right hand side of Eq. (61) is the diffusion 
coefficient of long thin rods, with D = D\\ = 2D±, including excluded volume 
corrections, with vo = 2l 2 /it. The second term on the right hand side of Eq. (61) 
arises from filament bundling driven at the rate a given in Eq. (26) and promotes 
density inhomogeneities. Equation (60) for the concentration couples to the Stokes 
equation, Eq. (57), with 

n*(c,n) = k B T a c(l + —\ + ra a a 5 ^^° c 2 , (62) 
V 7r / i62D 

and ^ 

°ii = ( 2r?0 + ~9W c ) w ^ ■ (63) 
In an isotropic active suspension there are no active contributions to the deviatoric 
part of the stress tensor, which has the form usual for passive rods [61]. There is, 
however, an active contribution to the pressure corresponding to the second term 
on the right hand side of Eq. (62). The first term on the right hand side of Eq. (62) 
is standard for passive rods. 

The homogeneous isotropic state in a quiescent suspension is characterized by 
v = and c = co- As discussed in the literature [30, 33, 36], the homogeneous state 
becomes unstable at high filament and motor concentration due to contractile effects 
generated by motor-induced filament bundling. Bundling is the main mechanism re- 
sponsible for the instability of both isotropic and ordered homogeneous states in 
quiescent suspensions. It is therefore instructive to display explicitly the details of 
this instability for the simple isotropic case. The examine the dynamics of fluctua- 
tions in the isotropic state we let c = Co + 8c and v = Sv in Eq. (60) and only keep 
terms of first order in the fluctuations. Incompressibility requires V • <5v = and the 
linearized equation for 5c is simply 

d t 8c = V(c )V 2 8c . (64) 

Expanding 8c in Fourier components, 8c(r, t) = = Ck(t)e lk r , one finds immedi- 
ately that the relaxation of the Fourier amplitudes, c k (t) = c k e^ Zc< ~ k ^ t , is controlled 
by a diffusive mode 

z c (k) = V(c )k 2 . (65) 

Density fluctuations become unstable when z c (k) < 0, corresponding to T>(co) < 
or c > cb, where 
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3D 3D 
Cb = -rz ^= — ~ -7-z — (66) 

is the concentration above which bundling overtakes diffusion. Using a ~ (b/l)uo, we 
can express the density cb in terms of the activity parameter fi a defined in Eq. (40) 
as cb = ^—(I'yp/a), where we have used D r = D/(6l 2 ). A possible location of this 
instability line in the phase diagram is shown in Fig. 4. 




Fig. 4. (color online) The phase diagram of homogeneous states for /i s = in the 
plane of filament density, Co, and motor activity /x a , as defined in Eq. (40), showing 
the location of the bundling instability at co = cb- The horizontal line at co = cjv 
for the isotropic-nematic transition crosses cip at fj, a fi x = 1/cjv- The cb line may lie 
above the cjvp — Cip line or cross through the iVand / states, as shown in the figure 
(l~fp/a = 0.1, ci4 = 50), depending on the value of l^p/a, a numerical parameter 
to leading order independent of ATP consumption rate. The instability of the / and 
N states is diffusive (dashed line), while the instability of the P state is oscillatory 
(dotted line). 



7.2 Nematic State 

The continuum variables describing the large-scale dynamics of an active nematic 
solution are the density and flow velocity of the solution and the concentration and 
alignment tensor of the filaments. For simplicity we consider only the case where 
there are no stationary apolar cross-linkers, i.e., m a — 0. In this case the transition 
from the isotropic to the nematic state occurs at the value cjv of passive suspensions. 
A finite fraction of stationary apolar crosslinkers lowers the value of the transition 
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density, as discussed in Section 6. In addition, it tends to stiffen all the liquid crystal 
elastic constants [36]. In the absence of external forces, the equations for filament 
concentration and alignment tensor are given by 

(d t + v ■ v) c = diDijdjC + diV Q dj {cQij) , (67) 



(d t + v v)(cQ, 3 ) = cF ij (k,Q)+H ij (c,Q) 



(68) 



The tensor Fij describes anisotropic convective and flow alignment effects and has 
the familiar form for passive nematic liquid crystals, as given in Eq. (14). At low 
density with the closure approximation described in Ref. [36] the alignment param- 
eter has the value A = l/(25o), with So the nematic order parameter defined in 
Eq. (10). The tensor Hij plays the role of the equilibrium molecular field for passive 
nematic liquid crystals, but it contains various active corrections. It is given by 



Hij(c, Q) ~ KV 2 (cQij) + K'[did k (cQj k ) + d 3 d k {cQ ik ) - Sijdkdi(cQ k i)] 



+dk(Ki jk idic) - 4D r (l - -^-^cQjj - D r a 4 c 3 Q k iQ k iQ 



(69) 



where 



3D 

2^(c,Q) = ^ 



1 + ( 1 - -S" ) — 

cat 



2 2\ 3c 4arh a c 



3D 



+ 
D 



fDv 4 . 



Kijki(c) 



cQij 
2arh a c 



D 2 

— (1 + v c) - -arh a c {5 ik 8ji + 8 jk S a 



■ SijS k i) 



am a c 
18 



(70) 
(71) 
(72) 
(73) 
(74) 



The last term on the right hand side of Eq. (69) , with 04 > has been introduced 
by hand. It arises from a quartic terms in the free energy of an equilibrium nematic 
and determines the magnitude of orientational order in passive rod solutions. It is 
apparent from the form of the various elastic constants in Eqs. (70-74) that bundling 
(described by the parameter a of Eq. (26)) always decreases the elastic constants of 
the nematic and therefore ultimately renders the uniform ordered state unstable. 

The flow velocity of the suspension is again obtained from Stokes' equation, 
Eq. (57), with a rods' contribution to the pressure given by 



n^(c,n) = kBT a c 



1+ 

+maQ 432^ C ( 5 ~ S ) ' 



+ kBT a ——u k iQ k i 
SbD 



(75) 



where the last term is new and arises from activity. The filament contribution to 
the deviatoric stress tensor is given by 
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= 2k B T a c(l - J1)q„ + rh a a^-c 2 ( > 
r rl 9 / 

+k B T a 



24D 



2 Ui i + 3 (uikQkj + ujkQki — SijUkiQklj^] ■ (76) 



Activity modifies the stress tensor of a nematic in two ways. The first term 
on the right hand side of Eq. (76) is equilibrium-like, in the sense that it can 
be obtained from the corresponding term in the stress tensor of passive rods, 
^passive _ 2ksT^l — -^^jQij by letting T — > T a (and replacing the transition 

density cjv by cin, when m s =^ 0). The second term on the right hand side of 
Eq. (76) is a truly nonequilibrium contribution. It was first proposed phenomeno- 
logically by Hatwalne and collaborators [50] who argued that an active element 
in solution behaves like a force dipole. Correlations among the axis of each dipole 
build up orientational order and yield active contributions to the stress tensor pro- 
portional to the orientational order parameter, Qij. Our microscopic derivation [57] 
yields an estimate for the coefficient of this term (undetermined, even in sign, in the 
phenomenological theory) and shows that the active cross-linkers yield contractile 
stresses (a > 0). Finally, the third term on the right hand side of Eq. (76) is the 
viscous contribution which has the standard form for a solution of rod-like filaments. 
Finally we note that active contributions proportional to the parameter (3 given in 
Eq. (27) do not appear in the hydrodynamics of the nematic phase. This is expected 
as terms proportional to /3 break the inversion symmetry of the ordered state and 
can only appear in a system with polar order. 



7.3 Polarized State 

The coarse-grained variables describing the dynamics of an active polarized suspen- 
sion are the density and flow velocity of the solution and the concentration and 
polarization of the filaments. As shown in Section 6, in a polarized state the align- 
ment tensor is slaved to the polarization field and it is not an independent continuum 
field. On the other hand, since our theory only considers terms that are quadratic in 
the fields, a nonzero value for |P| is only obtained by considering the coupled equa- 
tions for P to Qij and eliminating Qij in favor of P in the polarization equation 
to generate a term of order (P) 3 . To see, consider a filament density well into the 
polarized state, with c> cin and c > c/p, so that both coefficients a\ — 1 — c/c/p 
and (Z2 = 1 — c/cin in Eqs. (44) and (45) satisfy a\ < and d2 < 0. Setting the left 
hand side of Eqs. (44) and (45) to zero, we solve Eq. (45) for Qij to obtain 

QlJ = ^ { p i p i-\^P 2 ) ■ (77) 

This solution, substituted in Eq. (44), yields a term ~ P 2 Pi on the right hand side 
of Eq. (44) which has solution P 2 = (2a 1 a 2 )/(bib 2 c 2 ). 

The continuum equations for the polarized state are obtained by assuming that 
the alignment tensor relaxes on microscopic time scales to the form given by Eq. (77), 
which is then used to eliminate Qij in favor of P. With the exception of homoge- 
neous terms, such as the C((P) 3 ) term just described, this leads to a high density 
renormalization of the various coefficients in the continuum equations, but does not 
generate any new terms. For the sake of simplicity in the following we neglect this 
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renormalization and only keep those terms in the polarization equation generated 
by the coupling to the alignment tensor that have a qualitatively new structure. We 
also neglect all excluded volume corrections. The equation for filament concentration 
is given by 

ftc = -V ■ c(v - ^m a p C P) + di(V^{c)d jC ) 

-^arhadi^djiPiPj)] , (78) 

with 

23P.(c,P) = - afhatyij - arh^PP, + \SijP 2 ) ■ (79) 
The equation for the polarization vector has the form 

(ft + v ■ V) (cP) = i ( V x v) x (cP) + ^ [ Vv + ( V v) T ] • cP 

+H(c,P). (80) 

where H(c, P) generalizes the molecular field of equilibrium polar fluids [63] by 
including active contributions. It is given by 

H t (c, P) ~ - [D T - -yprhac + a 3 P 2 ] cPi + ^m a pcd t c - -£™a/?ft [c (PiPj - ^P 2 )] 
+ [d j K p d i (cP j )+d i K p d j (cP j j\ +d j K p d j (cP i ) 

-a J of Jfe ( c ,p)a fcC + 7 p^v 2 ( c p ! ) , (si) 



where 



f I \ D arha fool 

K P (c) = — — c, (82) 

K p( c ) = -S 7~ c > ( 83 ) 



^■*(c,P)=c 

The parameter 03 determines the value Po of the magnitude of the polarization in 
a quiescent (active) suspension, with P ( 2 = a'j,/[D r (c/cip — 1)]. 

In contrast to the case of the nematic, all three active mechanisms of motor- 
induced filament dynamics controlled by a, f3 and 7p appear in the hydrodynamic 
equations of polarized active suspension. Polarization sorting at a rate f3 ~ uo yields 
novel convective contributions in the first term on the right hand side of the equa- 
tion for the filament concentration, Eq. (78). In an equilibrium suspensions the 
filament concentration is convected with the suspension flow velocity, v. In an ac- 
tive polar suspension, in contrast, the filament concentration is convected with the 
effective velocity ~ v + m a /3P. The terms linear in the gradients proportional to 
f3 in the polarization equation are of similar origin. These terms were also incor- 
porated in the continuum description of self-propelled particles proposed by Simha 
and Ramaswamy. Bundling effects controlled by the rate a soften both the diffusion 
constant 23?. (c, P) in the concentration equation and the effective bend and splay 



Dvo am a 



)(PiS jk +S ij P k )+( 



nVvo , 2am a . ,. 

H — \Pjdik 



(84) 
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elastic constants K v and K' p of the polar fluid. Finally, the rotation rate 7p builds 
up polar order and controls the very existence of a polar state. 

For an incompressible suspension, the flow field v is obtained again from the 
Stokes equation, Eq. (57). The filament contribution to the pressure is given by 



i7 r P = k B T a c{ 1 + - 

TV 



k B T, 



U4D 



(85) 



The filament contribution to the deviatoric stress tensor of a polarized suspension 
is 



~r,P ~ ksT a 2 / o o 

= ma ° l2D C \ PiP * 2 



k B T a 
48Z> ' 



+ ■ 



k B T a c 2 b 2 
36D a 2 



UikPkPj + UjkPkPi — SijPktikiPi - UijP" 



. . r,k B T a 2 
+ma/3 4320 C 



diPj - 



(86) 



The first term on the right hand side of Eq. (86) is the active contribution to the 
stress tensor first discussed by Hatwalne and collaborators for a nematic suspen- 
sion [50]. The second and third term arise from the viscous coupling of filaments 
to the solvent. Finally the last term contains active contributions proportional to 
gradients of the polarization. These are controlled by the polarization sorting rate 
/3 ~ uo- Terms of these type are unique to the polarized state and vanish in a ne- 
matic suspension. They are expects to play an important role in renormalizing the 
rate of stress relaxation via reptation. 

Continuum equations for a polarized active suspension have been written down 
phenomenologically by several authors [49, 48, 51]. It is useful to make contact 
with this work. The phenomenological description can be recovered from our model 
by making a few simplifying approximations. An equation for the concentration of 
filaments of the form given in Eq. (78) was proposed by Ramaswamy and collab- 
orators [49, 51], although these authors neglected the diffusion terms, which play 
a crucial role in controlling the bundling instability of quiescent suspensions. The 
equation for the polarization vector P reduces to the form used by Voituriez et 
al. [48] and by Simha et al. [49, 51], if all terms containing higher order gradients of 
the concentration (PiV 2 c, Pi(Vc) 2 , P- Vftc, (P- Vc)(9»c)), as well as terms contain- 
ing both gradients of concentration and of polarization ((V ■ P)(9»c), (djc)(djPi), 
(djc)(diPj)) are neglected. With this approximation Eq. (80) becomes 



1 + Af 



{djVi)Pj 



1-Af 



(diVj)Pj 



-wic(P • V)Pi - w 2 cPi{V ■ P) + w 3 cd t \P\ 



+ 



^W4 + W 5 P 2 ^j 



W4+W 5 P ) Sij - W 6 Pi Pj 



djC 



+(K 1 -K 3 )d i V -P + K 3 V 2 P , 



(87) 



where r = ^pfh a c — D r > and P = r/as. Here the coefficients Wi have the 
form Wi = Cifkafi, with a numerical coefficients of order one. Note, however, that 
the terms proportional to Wi other than wi are equilibrium-like, in the sense that 
they could also be obtained from a polar contribution to the free energy of the form 
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. The wi term, in contrast, 



5F P = j r Bi5c(V-P)+B 2 c|P| 2 (V-P)+S 3 |P| 2 P-Vc+. 

is a true nonequilibrium contribution induced by activity and cannot be obtained 
from a free energy. All the remainder Wi's contain in general both equilibrium- 
like contributions determined by the Bi and nonequilibrium ones proportional to 
(3 ~ Ratp- Finally, K\ and K3 are the splay and bend elastic moduli, respectively, 
with 

Ki{c)= 7D_3n^c + n^c^ (gg) 

Ks(c) = f-^ + ^. (89) 

The first term on the right hand side of Eq. (87) guarantees the formation of a 
uniformly polarized state with |P| = Po- The next two terms are conventional cou- 
plings of liquid crystalline order and flow, with Ap the flow alignment parameter. 
Our low density calculation yields Ap = 1/2. The three terms on the second line 
are nonequilibrium terms analogous to those first written down by Toner and Tu in 
models of flocking [53, 54, 55]. The third line describes nonequilibrium changes in 
polarization driven by concentration gradients. Only the first of these terms (~ W4) 
is generally included in phenomenological theories. Equations (88) and (89) show 
that motor-induced filament bundling (~ a) softens both the splay and bend elastic 
constants, while polar crosslinkers (~ 7p) tend to stiffen them. Such effects, i.e. the 
dependence of elastic constants on the active elements are clearly beyond the scope 
of phenomenological theories with arbitrary elastic constants. In addition, the mi- 
croscopic derivation also provides contributions to the stress tensor which are higher 
order in gradients without the need for new unknown parameters, e.g. the expression 
for the stress tensor given in Eq. (86) contains novel contributions proportional to 
gradients of polarization that were not considered by other authors [49, 50, 46, 47] 
but whose microscopic origin is the same as those of lower order in the gradients. 
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A Appendix: Derivation of the rods' stress tensor 

We model very long, thin rods as rigid strings of spherical beads of diameter b <C I 
suspended in a fluid of viscosity rjo- We assume each rod consists of an odd number 
l/b = 2M + 1 of such beads, a sketched in Fig. 5. The beads on the a-th rod are 
indexed by an integer m that runs from —M to +M and the center of the m-th 
bead is at r a (m) = r a + mbu a , with r a = X)m rc "( m ) tne cen ter of mass of the 
rod. Momentum conservation at low Reynolds number is described by the Stokes 
equation 

M 

»7oV 2 v(r) - Vi7 - ^ *(r-r„(m))f£(m)=0, (90) 

ol m= — M 
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Fig. 5. (color online) The bead model of a rigid filament. 



where we have modeled each bead as a point-force on the fluid at the position of 
its center of mass. Faxen's theorem for the hydrodynamic force on a sphere in an 
inhomogeneous flow relates the force f^(m) exerted by the fluid on a bead to the 
flow velocity field vo(r) at the bead's position in the absence of that bead, according 
to 

f a (m) = — C&[v«(m) - v (r a (m))] , (91) 

where v a (m) = v a + mbuj a x u a is the velocity of the bead, with v a and u} a the 
center of mass and angular velocity of the rod, and £(, = Snrjob is the Stokes friction 
coefficient of a sphere of diameter b in an unbounded fluid of viscosity rjo- Using the 
linearity of the Stokes equation and the principle of superposition, the velocity of 
fluid at the position of the bead is given by 

v (r a (m)) = v(r a (m)) - H (r a (m) - r a (n)) ■ f£(n) , (92) 

where v(r) is the velocity of the fluid taking account of the presence of other rods 

and Hij(r) = (Su + rifj) is the Oseen tensor. Here the hydrodynamic in- 

teractions between beads on the same rod have been included explicitly, while the 
hydrodynamic coupling to other rods is implicitly taken into account in determining 
the flow velocity v(r). The force on bead m on the a— th rod is therefore given by 

fa(m) = -Cb [v Q (ra) - v (r a (m))l - ^ -. — - — :(S + u a ii a ) ■ ia{n) . (93) 

8 * — ' \n — m\ 

Now we take the limit I 3> b and introduce the continuous variable s — bm, where 
—1/2 < s < 1/2, so that r Q (s) — r a + su a . The hydrodynamic force per unit length, 
J 7 a( s )j satisfies the equation 

^i(*) = -C(vaW-v(r a ( a )))-| f ^—(S + u a u a )-T':(s'), (94) 

J \s — s' |>fe I I 

where £ s = 3ttt}o = Cb/b. Approximating 
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5{s - s) , (95) 



\s-s'\ \\s-s'\ 
where 

k[.[ H * s ' 1 ^ b) v^\ =2Hl/2b) ' (96) 

with 0(x) the Heaviside function, we obtain 

31n(f/2b) {s + ila{la) . jrh {s) „ _ Cs [Va(g) _ v (ra(s))] (9?) 

Since v a (s) = v a + sa; a x u a , then integrating equation (97) over s, we can obtain 
expressions for the hydrodynamic force, F^ = (J r ^(s)) s and torque = (u a s x 
JT^(s)) s at the center of mass of a rod, with (...) s = J^, 2 ds... as 

-C'^Ua) 1 1"! = v a - }<v(r a ( S ))) s , (98) 
-j-Ta = uf a - I~ 1 j(u a x sv(r a (s))) s , (99) 

where fo(u) = Cxfe - fiifl,-) + C||W C± = 2C|| = I^f), Cr = 3^6) and 
/ = Z 2 / 12. Performing a Taylor expansion of the fluid velocity about the center of 
mass, we obtain to lowest order in gradients 

-C _1 (ua) ■ Fa = v a - v(r a ) - I (iia ■ V) 2 v(r a ) + 0(V 4 ) , (100) 
~t* = u, a -u a x (u Q • V)v(r a ) + 0(V 3 ) . (101) 

Finally, we require that the hydrodynamic forces and torques be balanced by all 
other forces and torques on the rod. This gives 

F*=V a l/ M +fcflr Vlnc-fS, (102) 
tI = TZ a U ex + k B T a Tl a In c - , (103) 

where we have included contributions from fluctuations (non-equilibrium osmotic 
pressure), excluded volume interactions and active driving by the motors. Using 
Eqs. (97) and (100), we can calculate the hydrodynamic force per unit length on the 
rod as 

Fai{s) = C»j(Ua)[«aj ~ Vj{v a ) + S ((uf a X U a )j - (u ■ V a )Vj(r a )) 

-^(u.V a ) 2 ^(r a )] . (104) 

from which we obtain Eq. (32). Furthermore, defining the translational and rota- 
tional currents as 

JcM) = (^2v a 6(r-r a (t)^J (105) 
Jc(r,t) = (^2uf a 5(r-r a (t))j , (106) 
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the Smoluchowski equation (19) for the dynamics of c(r, u, t) follows. 

From equation (32), we perform the coarse-graining procedure to obtain the 
stress tensor. Retaining all terms of first order in gradients of the hydrodynamic 
fields, the pressure is 
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